Ejercicio 1

1.a) Se carga la informacion de los cajeros 108 al 116,

Caj108<-read.csv("cajero108.csv",header=F,dec=".",sep=";") 
Caj109<-read.csv("cajero109.csv",header=F,dec=".",sep=";") 
Caj110<-read.csv("cajero110.csv",header=F,dec=".",sep=";") 
Caj111<-read.csv("cajero111.csv",header=F,dec=".",sep=";") 
Caj112<-read.csv("cajero112.csv",header=F,dec=".",sep=";") 
Caj113<-read.csv("cajero113.csv",header=F,dec=".",sep=";") 
Caj114<-read.csv("cajero114.csv",header=F,dec=".",sep=";") 
Caj115<-read.csv("cajero115.csv",header=F,dec=".",sep=";") 
Caj116<-read.csv("cajero116.csv",header=F,dec=".",sep=";") 
dim(Caj108)

## [1]   1 122

dim(Caj109)

## [1]   1 243

dim(Caj110)

## [1]   1 396

dim(Caj111)

## [1]   1 251

dim(Caj112)

## [1]   1 185

dim(Caj113)

## [1] 215   1

dim(Caj114)

## [1] 215   1

dim(Caj115)

## [1] 215   1

dim(Caj116)

## [1] 215   1

Convertimos ahora cada una en una serie de tiempo, debemos tener cuidado puesto que las cinco primeras se encuentran como vectores filas,

Caj108<-t(Caj108)
Caj109<-t(Caj109)
Caj110<-t(Caj110)
Caj111<-t(Caj111)
Caj112<-t(Caj112)
Caj108<-ts(Caj108[,1],start=c(2008,1),freq=365)
Caj109<-ts(Caj109[,1],start=c(2008,1),freq=365)
Caj110<-ts(Caj110[,1],start=c(2008,1),freq=365)
Caj111<-ts(Caj111[,1],start=c(2008,1),freq=365)
Caj112<-ts(Caj112[,1],start=c(2008,1),freq=365)
Caj113<-ts(Caj113[,1],start=c(2008,1),freq=365)
Caj114<-ts(Caj114[,1],start=c(2008,1),freq=365)
Caj115<-ts(Caj115[,1],start=c(2008,1),freq=365)
Caj116<-ts(Caj116[,1],start=c(2008,1),freq=365)
str(Caj108)

##  Time-Series [1:122] from 2008 to 2008: 2440000 1248000 2499000 2079000 1335000 ...
##  - attr(*, "names")= chr [1:122] "V1" "V2" "V3" "V4" ...

str(Caj109)

##  Time-Series [1:243] from 2008 to 2009: 1989000 7780000 6859000 4802000 4733000 ...
##  - attr(*, "names")= chr [1:243] "V1" "V2" "V3" "V4" ...

str(Caj110)

##  Time-Series [1:396] from 2008 to 2009: 19222000 14495000 13983000 12565000 7124000 ...
##  - attr(*, "names")= chr [1:396] "V1" "V2" "V3" "V4" ...

str(Caj111)

##  Time-Series [1:251] from 2008 to 2009: 3562000 659000 2540000 1855000 1858000 ...
##  - attr(*, "names")= chr [1:251] "V1" "V2" "V3" "V4" ...

str(Caj112)

##  Time-Series [1:185] from 2008 to 2009: 4929000 1836000 1770000 1928000 4390000 ...
##  - attr(*, "names")= chr [1:185] "V1" "V2" "V3" "V4" ...

str(Caj113)

##  Time-Series [1:215] from 2008 to 2009: 24990 17780 18270 17260 20935 21810 9755 14620 10755 14000 ...

str(Caj114)

##  Time-Series [1:215] from 2008 to 2009: 1470 265 855 435 735 255 600 170 715 430 ...

str(Caj115)

##  Time-Series [1:215] from 2008 to 2009: 5450 3630 4710 3305 3590 4610 4630 2930 2470 3595 ...

str(Caj116)

##  Time-Series [1:215] from 2008 to 2009: 1405 1810 885 890 2020 1615 1395 1850 1145 990 ...

Ejercicio 2

Cargamos el paquete itsmr,

library(itsmr)
str(dowj)

##  num [1:78] 111 111 110 111 111 ...

str(strikes)

##  num [1:30] 4737 5117 5091 3468 4320 ...

str(Sunspots)

##  num [1:100] 101 82 66 35 31 7 20 92 154 125 ...

str(wine)

##  num [1:142] 464 675 703 887 1139 ...

Convertimos cada uno de los vectores de datos en series de tiempo. En este caso todos resultan ser vectores columnas y no exite informacion respecto a si las fechas se encuentran en listadas en orden creciente o decreciente por lo que se supondra que el orden de las fechas es exactamente el requerido

tdowj<-ts(dowj, start=c(1972, 241), freq=365)
plot(tdowj)

tstrikes<- ts(strikes, start=1951, end=1980, freq=1)
plot(tstrikes)

tSunspots<- ts(Sunspots, start=1870, end=1979, freq=1)
plot(tSunspots)

twine<- ts(wine, start=c(1980,1), end=c(1991,10), freq=12)
plot(twine)

str(tdowj)

##  Time-Series [1:78] from 1973 to 1973: 111 111 110 111 111 ...

str(tstrikes)

##  Time-Series [1:30] from 1951 to 1980: 4737 5117 5091 3468 4320 ...

str(tSunspots)

##  Time-Series [1:110] from 1870 to 1979: 101 82 66 35 31 7 20 92 154 125 ...

str(twine)

##  Time-Series [1:142] from 1980 to 1992: 464 675 703 887 1139 ...

Notemos que en el primer caso se utilizaron los par?metros \(1972\) y \(241\) para la fecha de iniciio de la serie, esto debido a que el 28 de agosto corresponde al dia 241 del a~no.

Ejercicio 3

Cargamos la tabla de datos DJTable.csv,

DJTable<-read.csv("DJTable.csv",header=T,dec=".",sep=";") 
tail(DJTable)

##              X    AA   AXP    BA   BAC   CAT  CSCO   CVX    DD   DIS    GE
## 247 11/01/2010 17.45 41.47 60.87 16.93 64.13 24.59 80.88 34.26 31.36 16.76
## 248 08/01/2010 17.02 41.95 61.60 16.78 60.34 24.66 79.47 33.94 31.88 16.60
## 249 07/01/2010 16.61 41.98 62.20 16.93 59.67 24.53 79.33 34.39 31.83 16.25
## 250 06/01/2010 16.97 41.49 59.78 16.39 59.43 24.42 79.63 34.04 31.82 15.45
## 251 05/01/2010 16.13 40.83 58.02 16.20 59.25 24.58 79.62 33.93 31.99 15.53
## 252 04/01/2010 16.65 40.92 56.18 15.69 58.55 24.69 79.06 34.26 32.07 15.45
##        HD   HPQ    IBM  INTC   JNJ   JPM    KO   MCD   MMM   MRK  MSFT
## 247 28.16 52.43 129.48 20.95 64.22 44.53 56.27 62.32 83.98 37.85 30.27
## 248 28.98 52.59 130.85 20.83 64.21 44.68 55.15 61.84 84.32 37.70 30.66
## 249 29.12 52.20 129.55 20.60 63.99 44.79 56.19 61.90 83.73 37.72 30.45
## 250 28.78 52.18 130.00 20.80 64.45 43.92 56.33 61.45 83.67 37.66 30.77
## 251 28.88 52.67 130.85 20.87 63.93 43.68 56.35 62.30 82.50 37.16 30.96
## 252 28.67 52.45 132.45 20.88 64.68 42.85 57.04 62.78 83.02 37.01 30.95
##       PFE    PG     T   TRV   UNH   UTX    VZ   WMT   XOM
## 247 18.83 60.20 26.97 48.54 32.92 72.16 31.88 54.21 70.30
## 248 18.68 60.44 27.10 48.56 32.70 70.63 31.75 53.33 69.52
## 249 18.53 60.52 27.30 48.63 33.01 70.49 31.73 53.60 69.80
## 250 18.60 60.85 27.61 47.94 31.79 70.19 31.92 53.57 70.02
## 251 18.66 61.14 28.44 48.63 31.48 70.56 33.34 53.69 69.42
## 252 18.93 61.12 28.58 49.81 31.53 71.63 33.28 54.23 69.15

Vemos que las fechas estan en orden decreciente por lo que debemos invertirlas, para ello,

head(DJTable)

##            X    AA   AXP    BA   BAC   CAT  CSCO   CVX    DD   DIS    GE
## 1 31/12/2010 15.39 42.92 65.26 13.34 93.66 20.23 91.25 49.88 37.51 18.29
## 2 30/12/2010 15.21 42.51 65.01 13.28 93.87 20.23 91.60 49.69 37.48 18.19
## 3 29/12/2010 15.13 42.86 65.05 13.31 93.78 20.25 91.37 50.02 37.60 18.27
## 4 28/12/2010 15.25 42.79 64.86 13.34 93.69 20.35 91.19 49.82 37.36 18.32
## 5 27/12/2010 15.23 43.05 64.75 13.27 94.07 20.16 90.12 49.63 37.48 18.19
## 6 23/12/2010 15.34 42.77 65.06 13.06 94.45 19.69 90.68 49.77 37.70 18.04
##      HD   HPQ    IBM  INTC   JNJ   JPM    KO   MCD   MMM   MRK  MSFT   PFE
## 1 35.06 42.10 146.76 21.03 61.85 42.42 65.77 76.76 86.30 36.04 27.91 17.51
## 2 34.86 42.26 146.67 21.02 61.94 42.23 65.50 76.76 86.54 36.01 27.85 17.49
## 3 34.89 42.32 146.52 20.94 62.13 42.36 65.45 76.99 86.76 36.21 27.97 17.60
## 4 35.09 42.25 145.71 20.88 62.05 42.61 65.36 76.43 86.74 36.20 28.01 17.59
## 5 35.24 41.82 145.34 20.84 61.93 42.67 65.07 76.43 87.01 36.23 28.07 17.49
## 6 35.09 41.74 145.89 20.84 62.25 42.08 65.58 76.96 86.47 36.29 28.30 17.61
##      PG     T   TRV   UNH   UTX    VZ   WMT   XOM
## 1 64.33 29.38 55.71 36.11 78.72 35.78 53.93 73.12
## 2 64.28 29.33 55.54 35.94 78.85 35.56 54.07 73.36
## 3 64.40 29.31 55.59 35.91 79.10 35.58 54.08 73.37
## 4 64.76 29.23 55.66 35.72 79.29 35.62 53.74 73.42
## 5 64.67 29.25 55.79 35.54 79.27 35.50 53.57 73.01
## 6 65.24 29.20 55.48 35.77 79.50 35.44 53.60 73.20

CSCO<-rev(DJTable$CSCO)
CVX<-rev(DJTable$CVX)
DD<-rev(DJTable$DD)

Por ultimo las convertimos en series de tiempo,

CSCO<-ts(CSCO, start=c(2010, 1), freq=365)
CVX<-ts(CVX, start=c(2010, 1), freq=365)
DD<-ts(DD, start=c(2010, 1), freq=365)
str(CSCO)

##  Time-Series [1:252] from 2010 to 2011: 24.7 24.6 24.4 24.5 24.7 ...

str(CVX)

##  Time-Series [1:252] from 2010 to 2011: 79.1 79.6 79.6 79.3 79.5 ...

str(DD)

##  Time-Series [1:252] from 2010 to 2011: 34.3 33.9 34 34.4 33.9 ...

Ejercicio 4

Cargamos la informacion del cajero 103,

Caj103<-read.csv("cajero103.csv",header=F,dec=".",sep=";")
dim(Caj103)

## [1] 1065    5

Debemos transformar a Caj103 en un vector columna, para ello

Caj103<-as.matrix(Caj103)  ### Pues Caj103 es un DataFrame
Caj103<-as.vector(Caj103)
length(Caj103)

## [1] 5325

Convertimos ahora el vector en una serie de tiempo

Caj103<-ts(Caj103,start=c(1998,1),freq=365) 
plot(Caj103)

Para verificar la normalidad graficamos el histograma para los residuos,

hist(diff(Caj103),prob=T,ylim=c(0,2.5e-07),col="red")
lines(density(diff(Caj103)),lwd=2) 
mu<-mean(diff(Caj103)) 
sigma<-sd(diff(Caj103)) 
x<-seq(-1e+07,1e+07,length=1000) 
y<-dnorm(x,mu,sigma) 
lines(x,y,lwd=2,col="blue")

Vemos de la grafica que el supuesto de normalidad se verifica aproximadamente bien.

Ejercicio 5

Realizamos el suavizamiento de la serie para \(a=4\), \(a=6\) y \(a=10\),

plot(Caj103, type="l")
st.1<- filter(Caj103,filter=rep(1/9,9)) 
st.2<-filter(Caj103,filter=rep(1/13,13)) 
st.3<-filter(Caj103,filter=rep(1/21,21)) 
lines(st.1,col="red") 
lines(st.2,col="purple") 
lines(st.3,col="blue")

Esta serie presenta mucha variabilidad, vemos como ella se reduce al hacer el suavizamiento o filtrado.

Ejercicio 6

Realizamos la dscomposicion de la serie utilizando la funcion stl. Para ello utilizamos la serie de tiempo Caj103 que ya ha sido generada,

DescCaj103<-stl(Caj103,s.window="periodic") 
head(DescCaj103$time.series)

##        seasonal   trend  remainder
## [1,] 1831668.00 3350290 -4553958.1
## [2,]  809017.04 3351466 -2998482.8
## [3,]  359299.40 3352641   179059.2
## [4,]   72981.77 3353817  -421798.8
## [5,]  -65002.53 3354993  -573990.2
## [6,] -571453.49 3356168  -163714.9

plot(DescCaj103)

Ejercicio 7

Cargamos el paquete itsmr y cargamos la tabla de datos deaths,

suppressMessages(library(itsmr))
str(deaths)

##  num [1:72] 9007 8106 8928 9137 10017 ...

Transformamos ahora la tabla en una serie de tiempo. Es importante recordar que los datos corresponden a muertes accidentales en Estados Unidos entre los a?os \(1973\) y \(1978\). Dado que la tabla tiene \(72\) entradas la frecuencia de losd datos es mensual, por lo tanto la frecuencia es igual a \(12\).

Deaths<-ts(deaths, start=c(1973,1), freq=12)

Ahora bien, calculamos primero el peridiograma

res<-spec.pgram(Deaths, log = "no")

res.ord<-order(res$spec,res$freq,decreasing=TRUE)
res.ord[1:4]

## [1]  6 12  1 30

Vemos que los m?ximos se encuentran en las posiciones \(6, 12 y 30\).

max1<-res$freq[6]
period1<-12/max1
max2<-res$freq[12]
period2<-12/max2
max3<-res$freq[30]
period3<-12/30
c(period1,period2,period3)

## [1] 12.0  6.0  0.4

res<-spec.pgram(Deaths, log = "no")
abline(v=max1, lty="dotted",col="red")
abline(v=max2, lty="dotted",col="blue")
abline(v=max3, lty="dotted",col="magenta")

Vemos que el primer periodo es \(12\), esto justifica el utilizar el modelo SARIMA con perdiodo \(12\).

Generamos el modelo y hacmeos la predicci?n. No quedaba claro del enunciado si una predicci?n correspond?a a predecir s?lo un nuevo dato o a realizar una predicci?n a voluntad. En este caso se eligi? predecir un a?o,

fit<-arima(Deaths,order=c(1,2,1),seasonal=list(order=c(2,1,2),period=12)) 
LH.pred<-predict(fit,n.ahead=12) 
layout(matrix(c(1,1,2,3),2,2,byrow=TRUE)) 
plot(Deaths,xlim=c(1973,1980),ylim=c(7000,12000),type="o") 
lines(LH.pred$pred,col="green",type="o") 
lines(LH.pred$pred+2*LH.pred$se,col="green",lty=3,type="o") 
lines(LH.pred$pred-2*LH.pred$se,col="green",lty=3,type="o")
acf(Deaths,main="Autocorrelaci?n?Simple",col="black",ylim=c(-1,1)) 
pacf(Deaths,main="Autocorrelaci?n?Parcial",col="black",ylim=c(-1,1))

Ejercicio 8

Cargamos la tabla de datos correspondiente al cajero \(101\). Recordemos que debemos trasponer este vector de datos para luego transformarlo en una serie de tiempo. La serie comienza el \(1\) de enero del a?o \(2008\), y tiene frecuencia diaria.

Caj101<-read.csv("Cajero101.csv",header=F,dec=".",sep=";") 
Caj101<-t(Caj101)#?Matriz?transpuesta 
Caj101<-ts(Caj101[,1],start=c(2008,1),freq=365)
plot(Caj101,type="o",col="blue")

Ahora bien, calculamos primero el peridiograma

res<-spec.pgram(Caj101, log = "no")

res.ord<-order(res$spec,res$freq,decreasing = TRUE)
res.ord[1:4]

## [1] 101 439 220 438

Vemos que los m?ximos se encuentran en las posiciones \(101, 439\) y \(220\).

max1<-res$freq[101]
period1<-365/max1
max2<-res$freq[439]
period2<-365/max2
max3<-res$freq[220]
period3<-365/max3
c(period1,period2,period3)

## [1] 15.207921  3.498861  6.981818

res<-spec.pgram(Caj101, log = "no")
abline(v=max1, lty="dotted",col="red")
abline(v=max2, lty="dotted",col="blue")
abline(v=max3, lty="dotted",col="magenta")

Vemos que el primer periodo es \(15\), este es el periodo m?s adecuado. Realizamos la predicci?n ajustando un modelo \(SARIMA(1,1,2)(2,1,2)\) con per?odo \(15\).

fit<-arima(Caj101,order=c(1,1,2),seasonal=list(order=c(2,1,2),period=15)) 
LH.pred<-predict(fit,n.ahead=365) 
layout(matrix(c(1,1,2,3), 2, 2, byrow = TRUE)) 
plot(Caj101,xlim=c(2008,2013),ylim=c(min(Caj101)-3500000,max(Caj101)+1000000),type="o") 
lines(LH.pred$pred,col="green",type="o") 
lines(LH.pred$pred+2*LH.pred$se,col="green",lty=3,type="o") 
lines(LH.pred$pred-2*LH.pred$se,col="green",lty=3,type="o")
acf(Caj101, main="Autocorrelaci?n Simple",col="black",ylim=c(-1,1)) 
pacf(Caj101,main="Autocorrelaci?n Parcial",col="black",ylim=c(-1,1))

Ahora bien, extraemos los ?ltimos dos meses de la serie,

Caj101<-Caj101[1461:1521]
Caj101<-ts(Caj101,start=c(2012,2),freq=365)
plot(Caj101,type="o",col="blue")

Ahora bien, calculamos el peridiograma

res<-spec.pgram(Caj101, log = "no")

res.ord<-order(res$spec,res$freq,decreasing = TRUE)
res.ord[1:4]

## [1]  2 18 19  1

Vemos que los m?ximos se encuentran en las posiciones \(2, 18\) y \(19\).

max1<-res$freq[2]
period1<-365/max1
max2<-res$freq[18]
period2<-365/max2
max3<-res$freq[19]
period3<-365/max3
c(period1,period2,period3)

## [1] 32.000000  3.555556  3.368421

res<-spec.pgram(Caj101, log = "no")
abline(v=max1, lty="dotted",col="red")
abline(v=max2, lty="dotted",col="blue")
abline(v=max3, lty="dotted",col="magenta")

Vemos que el primer periodo es \(32\), pero este periodo es demasiado grande dado que los datos corresponden s?lo a dos meses. Por esto se elegir? el siguiente periodo que es iguala \(3\).

fit<-arima(Caj101,order=c(7,1,7),seasonal=list(order=c(2,1,2),period=3))

## Warning in arima(Caj101, order = c(7, 1, 7), seasonal = list(order = c(2, :
## possible convergence problem: optim gave code = 1

LH.pred<-predict(fit,n.ahead=30) 
layout(matrix(c(1,1,2,3), 2, 2, byrow = TRUE)) 
plot(Caj101,xlim=c(2012,2012.25),ylim=c(min(Caj101)-5000000,max(Caj101)+10000000),type="o") 
lines(LH.pred$pred,col="green",type="o") 
lines(LH.pred$pred+2*LH.pred$se,col="green",lty=3,type="o") 
lines(LH.pred$pred-2*LH.pred$se,col="green",lty=3,type="o")
acf(Caj101, main="Autocorrelaci?n Simple",col="black",ylim=c(-1,1)) 
pacf(Caj101,main="Autocorrelaci?n Parcial",col="black",ylim=c(-1,1))

Vemos que en este caso la prediccion parece tener un comportamiento un mas dinamico y no solo a describir la media.

Ejercicio 9

Comenzaremos por generar las funciones que nos permitir?n calcular los diferentes tipos de error y calibrar el modelo de Holt-Winters.

La siguiente funci?n nos permitir? calcular diferentes tipos de error, Error Relativo, Error Cuadr?tico Medio, Total de Predicciones Sobre el Valor Real, Error Relativo para las Predicciones por Sobre el Valor Real.

La siguiente funci?n nos permitir? calibrar el modelo de Holt-Winters. Se modific? la funci?n presentada en el curso de modo de elegir el incremento como par?metro de la funci?n.

calibrar<-function(serie.aprendizaje,serie.testing, pot.incremento) {
  incremento<-10^{pot.incremento}
  error.c<-Inf
  alpha.i<-incremento  # alpha no puede ser cero
  while(alpha.i<=1) {
    beta.i<-0
    while(beta.i<=1) {
      gamma.i<-0
      while(gamma.i<=1) {
         mod.i<-HoltWinters(serie.aprendizaje,alpha=alpha.i,beta=beta.i,gamma=gamma.i)
         res.i<-predict(mod.i,n.ahead=length(serie.testing))
         error.i<-sqrt(ECM(res.i,serie.testing))
         if(error.i<error.c) {
           error.c<-error.i
           mod.c<-mod.i         
         }
         gamma.i<-gamma.i+incremento
      }
      beta.i<-beta.i+incremento
    }
    alpha.i<-alpha.i+incremento
  }  
  return(mod.c)
}

Cargamos el paquete itsmr y cargamos la tabla de datos deaths,

suppressMessages(library(itsmr))
str(deaths)

##  num [1:72] 9007 8106 8928 9137 10017 ...

Deaths<-ts(deaths, start=c(1973,1), freq=12)

Construimos las tablas de aprendizaje y prueba, donde esta ?ltima contiene los ?ltimos diez datos:

Deaths.train<-Deaths[1:62]
Deaths.train<-ts(Deaths.train, start=c(1973,1), freq=12)
Deaths.test<-Deaths[63:72]
Deaths.test<-ts(Deaths.test, start=c(1978,3), freq=12)

Calibramos el moldelo de Holt-Winters para la serie de datos de aprendizaje:

modelo.HW<-calibrar(Deaths.train,Deaths.test, -1)
modelo.HW

## Holt-Winters exponential smoothing with trend and additive seasonal component.
## 
## Call:
## HoltWinters(x = serie.aprendizaje, alpha = alpha.i, beta = beta.i,     gamma = gamma.i)
## 
## Smoothing parameters:
##  alpha: 0.1
##  beta : 0.3
##  gamma: 0.2
## 
## Coefficients:
##            [,1]
## a    8607.57551
## b      26.30563
## s1   -710.30464
## s2   -407.20322
## s3    216.10027
## s4    792.95226
## s5   1750.51793
## s6   1076.09047
## s7     55.89459
## s8    441.49007
## s9   -111.42316
## s10    47.73263
## s11  -876.26319
## s12 -1617.78781

Ahora bien, recordamos de la tarea anterior que el periodo asociado a estos datos es 12. Calibrando el modelo ARIMA mediante el uso de la funci?n \(auto.arima\) obtenemos:

suppressMessages(library(forecast))
auto.arima(Deaths.train)

## Series: Deaths.train 
## ARIMA(0,1,1)(0,1,1)[12]                    
## 
## Coefficients:
##           ma1     sma1
##       -0.4295  -0.4787
## s.e.   0.1368   0.2061
## 
## sigma^2 estimated as 116720:  log likelihood=-356.01
## AIC=718.03   AICc=718.56   BIC=723.7

Donde vemos que los par?metros ?ptimos a utilizar para ajustar el modelo ARIMA son \((0,1,1)(0,1,1)\), construyendo el modelo,

modelo.A<-arima(Deaths.train,order=c(0,1,1),seasonal=list(order=c(0,1,1),period=12))

Generamos \(10\) meses de predicciones utilizando ambos modelos. Primero, la predicci?n correspondiente al m?todo de Holt-Winters,

predict.HW<-predict(modelo.HW,n.ahead=10)
plot(predict.HW,xlim=c(1973,1979), ylim=c(6800, 11500), main = "Predicciones mediante modelo de Holt-Winters", col=2)
lines(Deaths.train,col=1)

y la predicci?n correspondiente al modelo SARIMA deperiodo \(12\) y par?metros \(()()\)

predict.A<-predict(modelo.A,n.ahead=10)
plot(predict.A$pred,xlim=c(1973,1979), ylim=c(6800, 11500), main = "Predicciones mediante modelo SARIMA", col = 2)
lines(Deaths.train,col=1)

Calculamos ahora el Error Cuadr?tico Medio asociado a cada predicci?n:

ECM.HW<-sqrt(ECM(predict.HW,Deaths.test))
ECM.A<-sqrt(ECM(predict.A$pred,Deaths.test))
Errores.ECM<-c(ECM.HW, ECM.A)
names(Errores.ECM)<-c("Holt-Winters", "SARIMA")
Errores.ECM

## Holt-Winters       SARIMA 
##     173.6803     494.5145

Vemos claramente que el modelo de Holt-Winters es mejor que el SARIMA, esto puesto que fue calibrado precisamente para minimizar el error cuadr?tico medio. Si utiliz?ramos otra medida de error podr?a darse una situaci?n diferente. De todas formas, la misma funci?n podria utilizarse para minimizar otros tipos de error.

Ahora bien, el Error Relativo para ambos modelos est? dado por:

ER.HW<-ER(predict.HW,Deaths.test)
ER.A<-ER(predict.A$pred,Deaths.test)
Errores.R<-c(ER.HW, ER.A)
names(Errores.R)<-c("Holt-Winters", "SARIMA")
Errores.R

## Holt-Winters       SARIMA 
##   0.01579028   0.04676109

En este caso observamos tambi?n que es el modelo de Holt-Winters el que tiene asociado un menor error. El mejor modelo basado en este criterio es entonces el modelo de Holt-Winters.

Ejercicio 10

Caj101<-read.csv("Cajero101.csv",header=F,dec=".",sep=";") 
Caj101<-t(Caj101)#?Matriz?transpuesta 
Caj101<-ts(Caj101[,1],start=c(2008,1),freq=365)
plot(Caj101,type="o",col="blue")

str(Caj101)

##  Time-Series [1:1521] from 2008 to 2012: 5409000 4703000 4243000 6395000 4645000 ...
##  - attr(*, "names")= chr [1:1521] "V1" "V2" "V3" "V4" ...

Repetimos los pasos seguidos en el ejercicio anterior, esto es, construimos las tablas de aprendizaje y prueba, donde esta ?ltima contiene los ?ltimos 31 d?as:

Caj101.train<-Caj101[1:1490]
Caj101.train<-ts(Caj101.train, start=c(2008,1), freq=365)
Caj101.test<-Caj101[1491:1521]
Caj101.test<-ts(Caj101.test, start=c(2012,30), freq=365)

Calibramos el moldelo de Holt-Winters para la serie de datos de aprendizaje:

modelo.HW<-calibrar(Caj101.train,Caj101.test, -1)
modelo.HW

## Holt-Winters exponential smoothing with trend and additive seasonal component.
## 
## Call:
## HoltWinters(x = serie.aprendizaje, alpha = alpha.i, beta = beta.i,     gamma = gamma.i)
## 
## Smoothing parameters:
##  alpha: 0.1
##  beta : 0
##  gamma: 0.5
## 
## Coefficients:
##               [,1]
## a     2271522.4035
## b         854.2779
## s1     -57574.8707
## s2    1530778.3887
## s3     314372.0384
## s4     350133.6202
## s5     434520.2718
## s6    -866898.5304
## s7   -1155317.7307
## s8    -484213.3089
## s9     280161.0677
## s10   -438313.8717
## s11     59972.5202
## s12    356671.6813
## s13  -1046435.2936
## s14   -744261.8615
## s15   1897087.0799
## s16   2310506.3525
## s17    198026.1161
## s18    -13928.9809
## s19    335633.4082
## s20   -732956.9313
## s21   -725233.2713
## s22   -410112.8102
## s23   -593490.9936
## s24  -1077907.6844
## s25      1175.2000
## s26    224114.4342
## s27   -522736.1947
## s28  -1194135.3761
## s29    -67387.9964
## s30   1186733.6711
## s31   1836692.9823
## s32   1030195.6703
## s33   1371681.5556
## s34   -252439.3857
## s35    -87414.9217
## s36    611178.6592
## s37   -195914.3419
## s38    -44589.1168
## s39    -63857.2025
## s40   -364342.5538
## s41  -1270778.2540
## s42  -1187134.2767
## s43   -685054.7330
## s44    698706.8585
## s45    383887.6042
## s46   1937734.5588
## s47    492766.5719
## s48    215962.4992
## s49    104531.3040
## s50    -77266.6927
## s51   -303137.5129
## s52   -579220.5384
## s53  -1239284.8362
## s54   -992544.0477
## s55  -1417499.8884
## s56  -1436511.5520
## s57   -823162.3830
## s58     20380.7885
## s59   3145156.1260
## s60    880246.2415
## s61   -473229.9273
## s62    919115.1425
## s63   -522243.8008
## s64   -672140.2023
## s65   -755060.4137
## s66   -471695.2599
## s67   -553080.1691
## s68     28628.6654
## s69  -1366367.6798
## s70   -885176.8804
## s71   -501991.8133
## s72   1326950.6865
## s73     89627.3050
## s74   1248168.1676
## s75    400711.5532
## s76   -843084.6366
## s77    640771.2888
## s78    137818.7581
## s79     64853.2447
## s80   -427491.5093
## s81    118539.4466
## s82   -431168.4284
## s83  -1537357.0308
## s84  -1275983.2206
## s85   -657962.6687
## s86    506166.1566
## s87  -1172892.1328
## s88  -1463692.8849
## s89    359786.0696
## s90   -610899.8290
## s91   -396678.7524
## s92    779821.2484
## s93   1155616.3219
## s94   -125192.6910
## s95    532472.1697
## s96    318850.0210
## s97  -1204803.0431
## s98  -1284592.0140
## s99   -818018.3294
## s100   760703.0856
## s101 -1650923.1224
## s102  -334419.0467
## s103 -1317132.8392
## s104 -2205939.0952
## s105  -214622.3190
## s106   905171.3820
## s107   440205.9355
## s108  -268751.6520
## s109  -861325.6783
## s110  -256006.2746
## s111 -1137266.4437
## s112 -1164405.8828
## s113 -1668804.4763
## s114 -1960150.9449
## s115 -1590926.7648
## s116 -1024037.9947
## s117   -32585.5337
## s118 -1596228.8421
## s119 -1218303.3608
## s120  1219534.8313
## s121  1050806.6307
## s122  1424433.0484
## s123  -988827.9473
## s124  -688577.0783
## s125  -895305.7415
## s126  -916996.6005
## s127  -281295.1088
## s128  -552943.6189
## s129  -576499.3522
## s130  -772528.9856
## s131   726592.9211
## s132 -1495977.2418
## s133  -562624.2538
## s134   552914.7338
## s135  1538591.6907
## s136  -159049.8362
## s137  -283167.8556
## s138  1339665.5038
## s139 -2061466.3671
## s140  -928108.0471
## s141  -523478.1660
## s142  -896653.6772
## s143  -980748.3795
## s144  -513257.5095
## s145  -339707.8401
## s146 -1744099.6188
## s147 -1519503.0303
## s148  -820309.5803
## s149  -263554.0469
## s150   537116.8739
## s151  1136530.7701
## s152   -51521.1827
## s153 -1364749.6578
## s154  -226667.5328
## s155    50031.6830
## s156  -493198.7546
## s157  -978270.1957
## s158  -101482.7305
## s159  -838761.1843
## s160 -1838976.7325
## s161 -1559949.4432
## s162  -908167.5582
## s163 -1401897.5196
## s164   379090.3464
## s165   928976.4844
## s166  1208967.0285
## s167  -945528.8514
## s168   459766.5331
## s169   608085.0186
## s170  -250707.6197
## s171  -723705.7659
## s172  -786469.0734
## s173  -634019.8944
## s174 -1221759.8176
## s175 -1483245.2442
## s176  -786543.2990
## s177  -984328.6882
## s178  -161510.2295
## s179  2019913.0496
## s180  1192017.7670
## s181  -297780.9912
## s182  -199565.5525
## s183   431219.5071
## s184   848933.3095
## s185   391177.1450
## s186   284714.4674
## s187  -130446.8790
## s188 -1200497.9256
## s189 -1247719.8615
## s190  -639552.0267
## s191  -739456.2777
## s192  -648597.7593
## s193   856233.7560
## s194  1165767.3946
## s195 -1026072.8610
## s196   143917.9949
## s197   100030.3470
## s198  1050055.5615
## s199   292348.0507
## s200   210470.9653
## s201  -402832.5320
## s202 -1054224.8482
## s203 -1190383.4195
## s204  -599996.1032
## s205 -1156013.5035
## s206 -1161336.1664
## s207   159195.0016
## s208  -349217.0999
## s209 -1551269.3683
## s210 -1170059.4173
## s211   525929.8008
## s212  1458934.4663
## s213  1724563.9160
## s214  -296335.1809
## s215   496886.7791
## s216  -746656.1970
## s217     2107.5115
## s218  -603547.2638
## s219  -385039.7857
## s220   -15107.2559
## s221    28646.1263
## s222  -554288.5892
## s223 -1435195.2973
## s224  -974264.4597
## s225   721813.4350
## s226  1690289.4023
## s227  1213113.6627
## s228   892992.6557
## s229   931711.4010
## s230  -721573.3761
## s231   200583.5348
## s232  -689677.6826
## s233  -214446.4390
## s234  -409734.4339
## s235  -446508.7964
## s236  -833243.4596
## s237 -1325624.7615
## s238 -1112237.0203
## s239 -1094700.6780
## s240   125554.9604
## s241  -362367.6669
## s242  1024495.6129
## s243   894839.1808
## s244 -1058108.4873
## s245  1006801.0201
## s246  1340146.5082
## s247   299666.3034
## s248    31225.3069
## s249   -51598.7067
## s250  -392685.8589
## s251 -1015582.6467
## s252 -1309427.7494
## s253  -813623.7402
## s254  -799871.0730
## s255  -391218.0587
## s256  1981980.1301
## s257    96744.7524
## s258 -1014427.6722
## s259  1599283.3060
## s260   844586.5270
## s261   248962.3226
## s262   -46900.6243
## s263  -258602.0729
## s264  -421730.2924
## s265  -951342.6386
## s266 -1042563.7235
## s267  -705933.0517
## s268  -761829.8545
## s269  -376244.7857
## s270   957613.0287
## s271   316849.6330
## s272  -878773.1371
## s273   680437.0065
## s274   645861.3716
## s275   995660.5598
## s276  1123135.7019
## s277   182811.5573
## s278   306185.4402
## s279 -1237698.3390
## s280  -863474.0847
## s281 -1122166.5140
## s282   -84173.3579
## s283   129642.1416
## s284  -364744.4754
## s285  -879511.8333
## s286 -1451715.1767
## s287  -156703.6883
## s288  1497951.1736
## s289   952692.4533
## s290   520701.5986
## s291   442103.6895
## s292 -1016546.4413
## s293 -1382460.3225
## s294  -734269.7635
## s295 -1163797.3848
## s296  -886136.7049
## s297  -826968.1253
## s298  -214610.7451
## s299  -447057.1887
## s300  -857719.3093
## s301  -423112.3321
## s302  -746496.8626
## s303   679345.0459
## s304  2794777.9016
## s305  2390695.0981
## s306  -149021.9089
## s307   118255.8809
## s308    26179.2061
## s309  -235481.6333
## s310  -375776.1404
## s311  -369359.2738
## s312  -426444.4102
## s313  -822409.5428
## s314 -1215928.3229
## s315 -1135169.6542
## s316  -631545.0012
## s317   453019.9995
## s318  1171214.2607
## s319  2753812.4202
## s320  -147868.5510
## s321   402593.8243
## s322   409111.0680
## s323  -232616.2614
## s324  -347333.8329
## s325   218011.7302
## s326   136338.5645
## s327  -311262.4753
## s328  -957009.3846
## s329  -345387.8165
## s330  -205125.5307
## s331   624859.5294
## s332  1579114.5696
## s333  2580632.4006
## s334  1747073.7404
## s335   785464.2236
## s336  2661759.3267
## s337  2669837.2039
## s338  1443485.1093
## s339  1278156.4800
## s340  2413080.7979
## s341   799938.1078
## s342  1288736.6400
## s343  1873821.8572
## s344  1162773.3790
## s345  1875684.0414
## s346  2325187.2224
## s347  3053383.1894
## s348  1998182.9231
## s349  2217575.3548
## s350  1668473.1205
## s351  2638687.4137
## s352  2669125.6905
## s353  1629548.1397
## s354  2632874.1944
## s355  1223346.7772
## s356  1236551.7183
## s357  1423710.8372
## s358  1548086.6157
## s359  -932974.6608
## s360   983174.3889
## s361   685792.3608
## s362   313177.5169
## s363  -642116.0335
## s364   413025.6429
## s365   258159.7706

Ahora bien, recordamos de la tarea anterior que el periodo asociado a estos datos es \(15\). Calibrando el modelo ARIMA mediante el uso de la funci?n \(auto.arima\) obtenemos:

suppressMessages(library(forecast))
auto.arima(Caj101.train)

## Series: Caj101.train 
## ARIMA(2,1,2)                    
## 
## Coefficients:
##           ar1     ar2      ma1      ma2
##       -0.4580  0.1426  -0.1807  -0.7834
## s.e.   0.0899  0.0483   0.0832   0.0823
## 
## sigma^2 estimated as 2.838e+12:  log likelihood=-23460.14
## AIC=46930.27   AICc=46930.31   BIC=46956.8

Donde vemos que los par?metros ?ptimos a utilizar para ajustar el modelo ARIMA son \((2,1,2)(2,1,2)\), construyendo el modelo,

modelo.A<-arima(Caj101.train,order=c(2,1,2),seasonal=list(order=c(2,1,2),period=15))

Generamos \(31\) d?as de predicciones utilizando ambos modelos. Primero, la predicci?n correspondiente al m?todo de Holt-Winters,

predict.HW<-predict(modelo.HW,n.ahead=31)
plot(Caj101.train,xlim=c(2008,2012.1), main = "Predicciones mediante modelo de Holt-Winters")
lines(predict.HW,col=2)

y la predicci?n correspondiente al modelo SARIMA de periodo \(15\) y par?metros \((2,1,2)(2,1,2)\)

predict.A<-predict(modelo.A,n.ahead=10)
plot(Caj101.train,xlim=c(2008,2012.1), main = "Predicciones mediante modelo SARIMA", col = 1)
lines(predict.A$pred,col=2)

Calculamos ahora el Error Cuadr?tico Medio asociado a cada predicci?n:

ECM.HW<-sqrt(ECM(predict.HW,Caj101.test))
ECM.A<-sqrt(ECM(predict.A$pred,Caj101.test))
Errores.ECM<-c(ECM.HW, ECM.A)
names(Errores.ECM)<-c("Holt-Winters", "SARIMA")
Errores.ECM

## Holt-Winters       SARIMA 
##    1501120.4     739513.6

Vemos claramente que el modelo SARIMA tiene un error cuadr?tico medio m?s bajo que el de Holt-Winters.

Ahora bien, el Error Relativo para ambos modelos est? dado por:

ER.HW<-ER(predict.HW,Caj101.test)
ER.A<-ER(predict.A$pred,Caj101.test)
Errores.R<-c(ER.HW, ER.A)
names(Errores.R)<-c("Holt-Winters", "SARIMA")
Errores.R

## Holt-Winters       SARIMA 
##    0.4070711    0.1319049

En este caso observamos que, al igual que en el caso del error cuadr?tico medio, el error relativo asociado al modelo SARIMA es menor que el asociado al m?todo de Holt-Winters

Basados en ambos errores, el modelo SARIMA parece ser el m?s adecuado para modelar estos datos. Se elije entonces este modelo para realizar la predicci?n para siete d?as. Debemos primero generar el modelo para todo el conjunto de datos,

modelo.A<-arima(Caj101,order=c(2,1,2),seasonal=list(order=c(2,1,2),period=15))

Ahora realizamos la predicci?n para siete d?as,

Caj101.pred<-predict(modelo.A,n.ahead=7) 
layout(matrix(c(1,1,2,3), 2, 2, byrow = TRUE)) 
plot(Caj101,xlim=c(2008,2012.1),ylim=c(0, max(Caj101)+1000000),type="o")
lines(Caj101.pred$pred,col=2)

Una idea, a priori, respecto a la calidad de esta predicci?n puede generarse gracias al c?lculo de los tipos de error asociados al modelo. Este c?lculo no utiliza los valores reales para los siete d?as predichos y por ende no es una medida necesariamente precisa. A posteriori, esto es, una vez conocido el valor de la serie para algunos de estos siete d?as, podr?a calcularse alg?n estimado del error para evaluar la predicci?n en base a esta nueva informaci?n y con ello decidir si esta predicci?n ha sido pertinente o no. Recordemos que lo ?ptimo seria actualizar las predicciones cada vez que contemos con nueva informaci?n puesto que mientras m?s datos sean predecidos, mayor es la varianza asociada a esta predicci?n y por ende mayor probabilidad de que el valor predicho se encuentre lejos del valor real.

Ejercicio 11

Cargamos la tabla de datos del Dow Jones (DJTable.csv), y la tabla traspuesta (DJTableTranspose.csv),

DJTable<-read.csv("DJTable.csv",header=T,dec=".",sep=";")
DJTableTranspose<-read.csv("DJTableTranspose.csv",header=T,dec=".",sep=";")

Realizamos un ACP para la tabla DJTable,

suppressMessages(library(FactoMineR))
res<-PCA(DJTable[,-1], scale.unit=TRUE, ncp=6, graph = FALSE)
plot(res, axes=c(1, 2), choix="ind", col.ind="Red",new.plot=TRUE)

plot(res, axes=c(1, 2), choix="var", col.var="blue",new.plot=TRUE)

y para la tabla traspuesta,

res<-PCA(DJTableTranspose[,-1], scale.unit=TRUE, ncp=6, graph = FALSE)
plot(res, axes=c(1, 2), choix="ind", col.ind="Red",new.plot=TRUE)

plot(res, axes=c(1, 2), choix="var", col.var="blue",new.plot=TRUE)

En el primer caso detectamos al menos tres grupos diferentes. Mirando el c?rculo de correlaciones vemos que cada uno corresponde a combinaciones de indicadores que tomaron valores altos ese d?a. Dado que en este caso los datos corresponden a valores de \(30\) indocadores econ?micos, no es sencillo realizar la interpretaci?n puesto que aparecen muchos efectos posiblemente mezclados.

En el segundo caso, se ha realizado un ACP tomando como atributos a la fecha en la que los correspondiente indicadores econ?micos tomaron los correspondientes valores. En este caso, pueden identificarse dos marcados grupos, el compuesto por el indicador \(IBM\) y otro compuesto por los dem?s. Analizando el c?rculo de correlaciones, concluimos que el primer grupo est? conformado por un indicador que toma valores altos en todos los d?as.

Ejercicio 12

Recordemos que los datos ya fueron cargados en memoria durante la resoluci?n del ejercico anterior. En este caso, realizamos clasificaci?n jer?rquica utilizando el paquete dtw,

suppressWarnings(suppressMessages(library(dtw)))
suppressWarnings(suppressMessages(library(rattle)))
clust.series <- hclust(dist(DJTable[,-1],method="DTW"), method="average")
par(mfrow=c(1,1))
plot(clust.series)

Graficamos ahora los centros de gravedad para los tres grupos,

centros<-centers.hclust(DJTable[,-1],clust.series,nclust=3,use.median=FALSE)
par(mfrow=c(3,1))
plot(centros[1,],type="o")
plot(centros[2,],type="o")
plot(centros[3,],type="o")

Las gr?ficas parecen ser bastante similares, dado que la escala de las gr?ficas es diferente, se mostrar?n en el mismo plot de modo de poder compararlas de forma m?s objetiva

par(mfrow=c(1,1))
plot(centros[1,],type="o")
lines(centros[2,],type="o", col=2)
lines(centros[3,],type="o", col=3)

Vemos que el primer cl?ster est? formado por aquellas fechas en las que los indicadores \(CAT\), \(CVX\), \(IBM\), \(KO\), \(MCD\) y \(MMM\) tomaron valores considerablemente m?s altos. El segundo cl?ster est? formado por aquellas fechas en las que los indicadores \(AA\), \(PFE\), \(UNH\), \(VZ\), \(WMT\) y \(XOM\) alcanzaron los valores m?s bajos, por ?ltimo, el tercer cl?ster est? formado por aquellas fechas en las que el indicador \(MSFT\) tom? el valor m?s alto y los indicadores \(DD\) y \(MCD\) el valor m?s bajo, los dem?s valores son bastante similares a los del segundo cl?ster.

Ejecricio 13

Aplicamos el m?todo de K-medias para las series de tiempo del Dow Jones con tres cl?sters,

grupos<-kmeans(DJTable[,-1],3,iter.max = 1000)
# grupos$cluster
# grupos$centers
par(mfrow=c(3,1))
plot(grupos$centers[1,],type="o")
plot(grupos$centers[2,],type="o")
plot(grupos$centers[3,],type="o")

Al igual que en el ejercicio anterior, las gr?ficas parecen ser bastante similares, por lo que se mostrar?n en el mismo plot de modo de poder compararlas de forma m?s objetiva

par(mfrow=c(1,1))
plot(centros[1,],type="o")
lines(centros[2,],type="o", col=2)
lines(centros[3,],type="o", col=3)

El an?lisis en este caso es exactamente igual al hecho en el ejercicio anterior ya que los grupos identificados al utlizar k-means fueron los mismos que utilizando clasificaci?n jer?rquica. Esto es, el primer cl?ster est? formado por aquellas fechas en las que los indicadores \(CAT\), \(CVX\), \(IBM\), \(KO\), \(MCD\) y \(MMM\) tomaron valores considerablemente m?s altos. El segundo cl?ster est? formado por aquellas fechas en las que los indicadores \(AA\), \(PFE\), \(UNH\), \(VZ\), \(WMT\) y \(XOM\) alcanzaron los valores m?s bajos, por ?ltimo, el tercer cl?ster est? formado por aquellas fechas en las que el indicador \(MSFT\) tom? el valor m?s alto y los indicadores \(DD\) y \(MCD\) el valor m?s bajo, los dem?s valores son bastante similares a los del segundo cl?ster.

Laboratorio - Series de tiempo

Sebastian Niklitschek-Soto

12 de mayo del 2017